JB-Scent-Net

#load network

#load network
setwd("~/Desktop/Scent-Net")
gs4_auth(email=TRUE)

   ℹ The googlesheets4 package is using a cached token for 'jbohey@uci.edu'.

net.long <- read_sheet("1plhxbVh5IQLNVMazqO4OvfM3X6vbll_Y6BTMJ94bkI4", na=c("","na", "NA")) %>%    
      mutate(pollinator_group= if_else(pollinator=="papilio_gothica", "butterfly", pollinator_group)) %>% #fix papilio_gothica from bee -> butterfly
      mutate(plant= recode(plant,senecio_tall_int = "senecio_integerrimus", agoceris_glauca = "agoseris_glauca", agoceris_aurantiaca="agoseris_aurantiaca", hydrophyllum_fendeleri="hydrophyllum_fendleri"))

   Auto-refreshing stale OAuth token.

   ✔ Reading from "caradonna_rmbl_interaction_networks_data_EDI".

   ✔ Range 'caradonna_rmbl_interaction_networks_data_EDI'.

net <- net.long %>% select(plant, pollinator, interactions) %>% pivot_wider(names_from = plant, values_from = interactions, values_fn= sum, values_fill= 0) %>%  #pivots data set longer so plants as columns and pollinators as rows, replaces NAs with zeros
    column_to_rownames("pollinator") #pollinators are in first column, makes them into rownames; #deletes the first pollinators column


net <- as.matrix(t(net)) # turns it from dataframe to matrix bc packages like matrix; t transposes (turns dataset sideways - rows -> columns; columns -> rows) 
polls <- colnames(net) #now pollinators are column names
plants<-rownames(net)
snet <- sortweb(net) #sort by row & column totals; sorted rows by rowsums and columns by column sums. Biggest numbers in top left of matrix

#load scent

#Go to comm_vols for code to make scent.net
#Loads in vol, filters with bouquet, takes mean of each voc for each species 
#From comm_vols.Rmd
scent.net <- read_csv("data/Volatiles_by_species_mean.csv") %>% column_to_rownames("species") %>% #species are now rownames
    as.matrix()   #network analysis likes matrix 

   Rows: 46 Columns: 356
   ── Column specification ────────────────────────────────────────────────────────
   Delimiter: ","
   chr   (1): species
   dbl (355): .alfa.-Copaene, .alpha.-Cubebene, .alpha.-Farnesene, .beta.-Bisab...
   
   ℹ Use `spec()` to retrieve the full column specification for this data.
   ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

scent.net <- scent.net[,colnames(scent.net)!="Cyclohexane, 1,3,5-trimethyl-2-octadecyl-"] #filtering this compound bc it has 0 interactions or emissions... 355 vs 354. breaks code later

#loads in chemical class and gets rid of ?
compound.class <-read_sheet("1BlDIOs4STe5ZTTPYE6qez1Nsey_ovEpe5m-SN8Q55gg", sheet= "chemical_class") %>% mutate(major_class= factor(major_class,levels=c("aliphatic", "monoterpene", "sesquiterpene", "nitrogenous", "sulfur", "benzenoid")) %>% fct_na_value_to_level("other")) %>% 
      filter(compound_long!="Cyclohexane, 1,3,5-trimethyl-2-octadecyl-") #filtering this compound bc it has 0 interactions or emissions... 355 vs 354

   ! Using an auto-discovered, cached token.
     To suppress this message, modify your code or options to clearly consent to
     the use of a cached token.
     See gargle's "Non-interactive auth" vignette for more details:
     <https://gargle.r-lib.org/articles/non-interactive-auth.html>
   ℹ The googlesheets4 package is using a cached token for 'jbohey@uci.edu'.
   ✔ Reading from "355 chemical compound class".
   ✔ Range ''chemical_class''.

#codes chemical class by colors
voc.class.colors <- setNames(brewer.pal(7,"Set3"), levels(compound.class$major_class))

##match paul’s plants with the plants we have vocs for

scentplants<-rownames(scent.net) #%>% tolower() %>%  str_replace(" ", "_")

paulsplants <- plants

#plants not in both datasets 
setdiff(scentplants, paulsplants) #gives us bonus plants

   [1] "aquilegia_caerulea" "frasera_speciosa"   "geum_aleppicum"

setdiff(paulsplants, scentplants) 

   [1] "hydrophyllum_fendleri" "phacelia_heterophylla"

bothplants <- intersect(scentplants, paulsplants) #plants found in amanda and paul data

Pollinator groups and family

polls.count <-count(net.long, pollinator_group,pollinator_family, pollinator) #counting number of rows for each pollinator
polls.class <- tibble(pollinator=polls) %>%  left_join(polls.count) %>% pull(pollinator_group)

   Joining with `by = join_by(pollinator)`

polls.class.dataframe <- tibble(pollinator=polls) %>%  left_join(polls.count)

   Joining with `by = join_by(pollinator)`

polls.family <- tibble(pollinator=polls) %>%  left_join(polls.count) %>% pull(pollinator_family)

   Joining with `by = join_by(pollinator)`

polls.order.all <- polls.class

#net.grp.all has all of pauls 45 plant species 
#net.grp only has plant species found in Amanda's and paul's datasets
net.grp.all <- aggregate(.~ polls.class, as.data.frame(t(net)), sum) #lump network by pollinator class; plants columns, pollinators rows
rownames(net.grp.all) <- net.grp.all[,1] #make plants rows, pollinators column names
net.grp.all[,1] <- NULL
net.grp.all <- as.matrix(t(net.grp.all)) #transposes matrix
net.grp<- net.grp.all[bothplants,] 

net.cut <- net[bothplants,]

Plant-pollinator network

##Pollinator groups

#Tells you plants major pollinator
plants.majorpoll <- setNames(factor(colnames(net.grp)[apply(net.grp, 1, which.max)]), rownames(net.grp))

#sets colors for ALL 5 pollinator groups
pcols.all <- setNames(brewer.pal(5,"Set2"), colnames(net.grp))

#Four colors for each major pollinator group
pcols <- pcols.all[ levels(plants.majorpoll)]

Proportion of visits by pollinator group

#gives % bee pollination for each plant species 
plants.bees <- net.grp[,"bee"]/rowSums(net.grp)
plot(sort(plants.bees))

scent.net.cut <- scent.net[rownames(scent.net)%in% paulsplants, ] #cuts scent to match what paul has pollinator data for

# Cap of plants.bees
plant.bee.cap <- capscale(sqrt(scent.net.cut)~plants.bees)
anova(plant.bee.cap)

   Permutation test for capscale under reduced model
   Permutation: free
   Number of permutations: 999
   
   Model: capscale(formula = sqrt(scent.net.cut) ~ plants.bees)
            Df Variance      F Pr(>F)
   Model     1   505636 1.0525  0.355
   Residual 41 19696248

#p=0.348
plot(plant.bee.cap, type="text")

#gives % fly pollination for each plant species 
plants.fly <- net.grp[,"fly"]/rowSums(net.grp)
plot(sort(plants.fly))

sort(plants.fly)

          agoseris_aurantiaca      amelanchier_alnifolia 
                  0.000000000                0.000000000 
         castilleja_sulphurea            gentiana_parryi 
                  0.000000000                0.000000000 
          ipomopsis_aggregata       lathyrus_lanszwertii 
                  0.000000000                0.000000000 
               mahonia_repens          mertensia_ciliata 
                  0.000000000                0.000000000 
         ranunculus_inamoenus               rosa_woodsii 
                  0.000000000                0.000000000 
              vicia_americana            viola_praemorsa 
                  0.000000000                0.000000000 
      delphinium_nuttallianum     hydrophyllum_capitatum 
                  0.001774623                0.010204082 
       campanula_rotundifolia       mertensia_fusiformis 
                  0.011428571                0.011627907 
     erythronium_grandiflorum      heliomeris_multiflora 
                  0.019230769                0.033676658 
            sedum_lanceolatum        heterotheca_villosa 
                  0.035398230                0.044919786 
           lomatium_dissectum       taraxacum_officinale 
                  0.076923077                0.102841678 
                senecio_serra       senecio_integerrimus 
                  0.125000000                0.174044876 
            arenaria_congesta           boechera_stricta 
                  0.178947368                0.181818182 
           erigeron_speciosus      solidago_multiradiata 
                  0.202226345                0.205357143 
         claytonia_lanceolata helianthella_quinquenervis 
                  0.244186047                0.267904509 
          fragaria_virginiana              linum_lewisii 
                  0.333333333                0.333333333 
          potentilla_hippiana       erigeron_flagellaris 
                  0.397425583                0.432900433 
          potentilla_gracilis       eriogonum_umbellatum 
                  0.433962264                0.442477876 
    pseudocymopterus_montanus                draba_aurea 
                  0.490196078                0.564285714 
         eriogonum_subalpinum  androsace_septentrionalis 
                  0.575000000                0.630000000 
               galium_boreale       achillea_millefolium 
                  0.875000000                0.942622951 
              agoseris_glauca 
                  1.000000000

#gives % butterfly pollination for each plant species 
plants.butterfly <- net.grp[,"butterfly"]/rowSums(net.grp)
plot(sort(plants.butterfly))

sort(plants.butterfly)

              agoseris_glauca      amelanchier_alnifolia 
                  0.000000000                0.000000000 
    androsace_septentrionalis     campanula_rotundifolia 
                  0.000000000                0.000000000 
         castilleja_sulphurea                draba_aurea 
                  0.000000000                0.000000000 
     erythronium_grandiflorum        fragaria_virginiana 
                  0.000000000                0.000000000 
              gentiana_parryi     hydrophyllum_capitatum 
                  0.000000000                0.000000000 
                linum_lewisii         lomatium_dissectum 
                  0.000000000                0.000000000 
               mahonia_repens          mertensia_ciliata 
                  0.000000000                0.000000000 
         mertensia_fusiformis  pseudocymopterus_montanus 
                  0.000000000                0.000000000 
         ranunculus_inamoenus               rosa_woodsii 
                  0.000000000                0.000000000 
            sedum_lanceolatum              senecio_serra 
                  0.000000000                0.000000000 
              vicia_americana            viola_praemorsa 
                  0.000000000                0.000000000 
        heliomeris_multiflora        potentilla_hippiana 
                  0.003934730                0.006436042 
      delphinium_nuttallianum       taraxacum_officinale 
                  0.006654836                0.008119080 
         achillea_millefolium helianthella_quinquenervis 
                  0.012295082                0.013262599 
        solidago_multiradiata        potentilla_gracilis 
                  0.013392857                0.018867925 
          heterotheca_villosa       claytonia_lanceolata 
                  0.022887701                0.023255814 
         lathyrus_lanszwertii          arenaria_congesta 
                  0.025000000                0.057142857 
         senecio_integerrimus       erigeron_flagellaris 
                  0.065494239                0.089466089 
         eriogonum_subalpinum             galium_boreale 
                  0.125000000                0.125000000 
         eriogonum_umbellatum         erigeron_speciosus 
                  0.150442478                0.179962894 
          ipomopsis_aggregata           boechera_stricta 
                  0.217391304                0.363636364 
          agoseris_aurantiaca 
                  1.000000000

#gives % hummingbird pollination for each plant species 
plants.hummingbird <- net.grp[,"hummingbird"]/rowSums(net.grp)
plot(sort(plants.hummingbird))

sort(plants.hummingbird)

         achillea_millefolium        agoseris_aurantiaca 
                  0.000000000                0.000000000 
              agoseris_glauca      amelanchier_alnifolia 
                  0.000000000                0.000000000 
    androsace_septentrionalis          arenaria_congesta 
                  0.000000000                0.000000000 
       campanula_rotundifolia       castilleja_sulphurea 
                  0.000000000                0.000000000 
         claytonia_lanceolata                draba_aurea 
                  0.000000000                0.000000000 
         erigeron_flagellaris         erigeron_speciosus 
                  0.000000000                0.000000000 
         eriogonum_subalpinum       eriogonum_umbellatum 
                  0.000000000                0.000000000 
     erythronium_grandiflorum        fragaria_virginiana 
                  0.000000000                0.000000000 
               galium_boreale            gentiana_parryi 
                  0.000000000                0.000000000 
   helianthella_quinquenervis      heliomeris_multiflora 
                  0.000000000                0.000000000 
          heterotheca_villosa     hydrophyllum_capitatum 
                  0.000000000                0.000000000 
         lathyrus_lanszwertii              linum_lewisii 
                  0.000000000                0.000000000 
           lomatium_dissectum             mahonia_repens 
                  0.000000000                0.000000000 
            mertensia_ciliata       mertensia_fusiformis 
                  0.000000000                0.000000000 
          potentilla_gracilis        potentilla_hippiana 
                  0.000000000                0.000000000 
    pseudocymopterus_montanus       ranunculus_inamoenus 
                  0.000000000                0.000000000 
                 rosa_woodsii          sedum_lanceolatum 
                  0.000000000                0.000000000 
                senecio_serra      solidago_multiradiata 
                  0.000000000                0.000000000 
         taraxacum_officinale            vicia_americana 
                  0.000000000                0.000000000 
              viola_praemorsa       senecio_integerrimus 
                  0.000000000                0.002425713 
             boechera_stricta    delphinium_nuttallianum 
                  0.090909091                0.224933452 
          ipomopsis_aggregata 
                  0.608695652

#gives % wasp pollination for each plant species 
plants.wasp <- net.grp[,"wasp"]/rowSums(net.grp)
plot(sort(plants.wasp))

sort(plants.wasp)

         achillea_millefolium        agoseris_aurantiaca 
                 0.0000000000               0.0000000000 
              agoseris_glauca      amelanchier_alnifolia 
                 0.0000000000               0.0000000000 
    androsace_septentrionalis          arenaria_congesta 
                 0.0000000000               0.0000000000 
             boechera_stricta     campanula_rotundifolia 
                 0.0000000000               0.0000000000 
         castilleja_sulphurea       claytonia_lanceolata 
                 0.0000000000               0.0000000000 
      delphinium_nuttallianum                draba_aurea 
                 0.0000000000               0.0000000000 
         erigeron_flagellaris         erigeron_speciosus 
                 0.0000000000               0.0000000000 
         eriogonum_umbellatum   erythronium_grandiflorum 
                 0.0000000000               0.0000000000 
          fragaria_virginiana             galium_boreale 
                 0.0000000000               0.0000000000 
              gentiana_parryi helianthella_quinquenervis 
                 0.0000000000               0.0000000000 
        heliomeris_multiflora     hydrophyllum_capitatum 
                 0.0000000000               0.0000000000 
          ipomopsis_aggregata       lathyrus_lanszwertii 
                 0.0000000000               0.0000000000 
                linum_lewisii         lomatium_dissectum 
                 0.0000000000               0.0000000000 
               mahonia_repens          mertensia_ciliata 
                 0.0000000000               0.0000000000 
         mertensia_fusiformis        potentilla_hippiana 
                 0.0000000000               0.0000000000 
    pseudocymopterus_montanus       ranunculus_inamoenus 
                 0.0000000000               0.0000000000 
                 rosa_woodsii          sedum_lanceolatum 
                 0.0000000000               0.0000000000 
         senecio_integerrimus              senecio_serra 
                 0.0000000000               0.0000000000 
         taraxacum_officinale            vicia_americana 
                 0.0000000000               0.0000000000 
              viola_praemorsa        heterotheca_villosa 
                 0.0000000000               0.0005347594 
        solidago_multiradiata        potentilla_gracilis 
                 0.0178571429               0.0377358491 
         eriogonum_subalpinum 
                 0.0750000000

##Heatmap1

#pp_heatmap
heatmaply(net, scale="none", showticklabels = c(F,T), margins= c(0,NA,0,0),
          scale_fill_gradient_fun = ggplot2::scale_fill_gradient(low = "white", high = "#000050", trans="log1p"),
          distfun=vegdist, 
          hclustfun=function(x) hclust(x, method="mcquitty"), 
          row_side_colors = data.frame(P=plants.majorpoll[rowSums(net)>0]),
          row_side_palette=function(n) pcols, 
          col_side_colors=data.frame(PollinatorGroup=polls.class), 
          col_side_palette=function(n) pcols.all)

   Warning in doTryCatch(return(expr), name, parentenv, handler): unable to load shared object '/Library/Frameworks/R.framework/Resources/modules//R_X11.so':
     dlopen(/Library/Frameworks/R.framework/Resources/modules//R_X11.so, 0x0006): Library not loaded: /opt/X11/lib/libSM.6.dylib
     Referenced from: <FFA47D77-8F35-36FC-B0E5-38351B8D9512> /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/modules/R_X11.so
     Reason: tried: '/opt/X11/lib/libSM.6.dylib' (no such file), '/System/Volumes/Preboot/Cryptexes/OS/opt/X11/lib/libSM.6.dylib' (no such file), '/opt/X11/lib/libSM.6.dylib' (no such file), '/Library/Frameworks/R.framework/Resources/lib/libSM.6.dylib' (no such file), '/Library/Java/JavaVirtualMachines/jdk-11.0.18+10/Contents/Home/lib/server/libSM.6.dylib' (no such file)

##Matrix plot

#interaction matrix plot
visweb(net, type="nested", circles=T, boxes=F, circle.max=1.2)

##Web plot

#fix names of species (plants rows, polls cols)
net.nice.labels <- net
rownames(net.nice.labels) <- str_to_sentence(str_replace_all(rownames(net), "_"," "))
colnames(net.nice.labels) <- str_to_sentence(str_replace_all(colnames(net), "_"," "))

#interaction web plot
par(mar=c(0,0,0,0))
plant.pollinator.network <- plotweb(sqrt(net.nice.labels), empty=F, text.rot=90, method="cca", col.high=pcols.all[polls.order.all], col.low=pcols[plants.majorpoll], col.interaction=pcols.all[polls.order.all], bor.col.interaction=pcols.all[polls.order.all], bor.col.high=pcols.all[polls.order.all], bor.col.low=pcols[plants.majorpoll])

library(ggplot2)
ggsave(plant.pollinator.network, file='plant.pollinator.networkplot.jpg', width=130, height=105, units= 'mm', dpi=300)

#Network topography

#computes network statistics (like nestedness)
#networkstats <- networklevel(net)
#write_csv(enframe(networkstats), "data/p-p network stats.csv")

#computes network statistics for voc-pollinator network
#scentnetworkstats<- networklevel(scentpollnet)
#write_csv(enframe(scentnetworkstats), "data/scent network stats.csv")
scentnet.stats <- read_csv("data/scent network stats.csv")

   Rows: 48 Columns: 2
   ── Column specification ────────────────────────────────────────────────────────
   Delimiter: ","
   chr (1): name
   dbl (1): value
   
   ℹ Use `spec()` to retrieve the full column specification for this data.
   ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

#the network statistics null model: TODO- broken
#network.nullmodel<-nullmodel(scentpollnet) 
#nullmodel.stats<- map(network.nullmodel, ~networklevel(.x )) #network.nullmodel[1:2]; index=c("NODF", "weighted NODF", "connectance", "H2", "modularity Q")

#nullmodel.stats %>% map_dfr(~enframe(.x), .id = "run") %>% write_csv("data/scent network null model statistics.csv")
nullmodel.stats <- read_csv("data/scent network null model statistics.csv")

   Rows: 48000 Columns: 3
   ── Column specification ────────────────────────────────────────────────────────
   Delimiter: ","
   chr (1): name
   dbl (2): run, value
   
   ℹ Use `spec()` to retrieve the full column specification for this data.
   ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

#idea - pivot wider so network properties are the columns, rows are run number, then average by values in column -> pivot longer or store in new variable
#nullmodel.stats %>% pivot_wider(names_from = "Name", values_from = "Value") %>% column_to_rownames("stats") %>% as.matrix()

null.NODF <- read_csv("data/scent network null model statistics.csv") %>% filter(name=="NODF")

   Rows: 48000 Columns: 3
   ── Column specification ────────────────────────────────────────────────────────
   Delimiter: ","
   chr (1): name
   dbl (2): run, value
   
   ℹ Use `spec()` to retrieve the full column specification for this data.
   ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

ggplot(null.NODF, aes(x=value))+ geom_histogram() + geom_vline(aes(xintercept=filter(scentnet.stats, name=="NODF") %>% pull(value)))

   `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Scent network

##Heatmap2

scent.net.cut <- scent.net[rownames(scent.net)%in% paulsplants, ] #cuts scent to match what paul has pollinator data for
plants.majorpoll.cut <- plants.majorpoll[rownames(scent.net.cut)] #cuts Pauls data to match Amanda's

plants.bees.cut<- plants.bees[rownames(scent.net.cut)]

heatmaply(sqrt(scent.net.cut), scale="none", showticklabels = c(F,T), margins= c(0,NA,0,0),
          scale_fill_gradient_fun = ggplot2::scale_fill_gradient(low = "white", high = "#000050", trans="log1p"),
          distfun=vegdist, 
          hclustfun=function(x) hclust(x, method="mcquitty"), 
          row_side_colors = data.frame(major.polls=plants.majorpoll.cut),
          row_side_palette=function(n) pcols )

   Warning in distfun(x): you have empty rows: their dissimilarities may be
                    meaningless in method "bray"

          #col_side_colors=data.frame(PollinatorGroup=scent.net), #put compound class here
        #  col_side_palette=function(n) pcols.all)

#Distance Matrix (for bray curtis distance) #Mantel test

#scent-plant network 
distance.matrix <- vegdist(scent.net.cut)
View(as.matrix(distance.matrix)) #tells you how different the smells are. If 0 there is no difference, scents are the same; 1 = no compounds in common

#plant-pollinator network
#how similar are the plants in terms of their pollinator array/composition
distance.matrix.pp <- vegdist(net.cut)
View(as.matrix(distance.matrix.pp))

#Mantel Test
#can do correlation, if scent similarity predicts pollinator --> mantel test
mantel(distance.matrix,distance.matrix.pp)

   
   Mantel statistic based on Pearson's product-moment correlation 
   
   Call:
   mantel(xdis = distance.matrix, ydis = distance.matrix.pp) 
   
   Mantel statistic r: 0.02341 
         Significance: 0.297 
   
   Upper quantiles of permutations (null model):
      90%    95%  97.5%    99% 
   0.0636 0.0844 0.1041 0.1341 
   Permutation: free
   Number of permutations: 999

#TODO: group by family then rerun vegdist

distance.matrix.pollgroup <- vegdist(net.grp)
mantel(distance.matrix, distance.matrix.pollgroup) 

   
   Mantel statistic based on Pearson's product-moment correlation 
   
   Call:
   mantel(xdis = distance.matrix, ydis = distance.matrix.pollgroup) 
   
   Mantel statistic r: -0.06377 
         Significance: 0.899 
   
   Upper quantiles of permutations (null model):
      90%    95%  97.5%    99% 
   0.0746 0.0954 0.1154 0.1374 
   Permutation: free
   Number of permutations: 999

#can't predict which pollinator group will visit plant based on VOCS


#relative amounts of visits from pollinator groups 
distance.matrix.pollgroup.relative <- vegdist(decostand(net.grp, method="tot"))
mantel(distance.matrix, distance.matrix.pollgroup.relative) 

   
   Mantel statistic based on Pearson's product-moment correlation 
   
   Call:
   mantel(xdis = distance.matrix, ydis = distance.matrix.pollgroup.relative) 
   
   Mantel statistic r: -0.1527 
         Significance: 0.985 
   
   Upper quantiles of permutations (null model):
     90%   95% 97.5%   99% 
   0.107 0.140 0.170 0.217 
   Permutation: free
   Number of permutations: 999

#similarity of the relative amounts of scent of each plant
distance.matrix.relative <- vegdist(decostand(scent.net.cut, method="tot"))

#is similarity of the relative amounts of scent of each plant correlated with the similarity in relative visits from each pollinator group
mantel(distance.matrix.relative, distance.matrix.pollgroup.relative)

   
   Mantel statistic based on Pearson's product-moment correlation 
   
   Call:
   mantel(xdis = distance.matrix.relative, ydis = distance.matrix.pollgroup.relative) 
   
   Mantel statistic r: -0.1544 
         Significance: 0.985 
   
   Upper quantiles of permutations (null model):
     90%   95% 97.5%   99% 
   0.110 0.143 0.172 0.205 
   Permutation: free
   Number of permutations: 999

#No correlation between scent similarity and similarity of the proportion of visits by pollinators or pollinator group (p> 0.1)


#not grouped by pollinator class but still relative
distance.matrix.pp.relative <- vegdist(decostand(net.cut, method="tot"))
mantel(distance.matrix.relative, distance.matrix.pp.relative)

   
   Mantel statistic based on Pearson's product-moment correlation 
   
   Call:
   mantel(xdis = distance.matrix.relative, ydis = distance.matrix.pp.relative) 
   
   Mantel statistic r: 0.05562 
         Significance: 0.192 
   
   Upper quantiles of permutations (null model):
     90%   95% 97.5%   99% 
   0.078 0.104 0.118 0.139 
   Permutation: free
   Number of permutations: 999

##Web plot VOCs-Plants

#interaction web plot
#plant-voc bipartite
plotweb(sqrt(scent.net), text.rot=90, method="cca", col.high=voc.class.colors[compound.class$major_class], bor.col.high=voc.class.colors[compound.class$major_class], col.low=pcols[plants.majorpoll], bor.col.low=pcols[plants.majorpoll], col.interaction=pal[5], bor.col.interaction=pal[5])

#Plant-scent-pollinator network ##NMDS

scent.net.cut <- scent.net[rownames(scent.net)%in% paulsplants, ]
scent.mds <- metaMDS(sqrt(scent.net.cut), autotransform = F, trace=F)
scent.op <- ordiplot(scent.mds, type="n")
points(scent.op, what="species", pch=3, col="grey50")
points(scent.op, what="sites", cex=1.5, pch=19, col=pcols[plants.majorpoll])
text(scent.op, what="sites", cex=0.8, col=pcols[plants.majorpoll])

#are scents of plants with different major pollinators different? no... 
#do phylogeny, plants in same genus will cluster together need to account for phylogeny

#cap

scent.cap <- capscale(sqrt(scent.net.cut)~plants.majorpoll.cut)
anova(scent.cap)

   Permutation test for capscale under reduced model
   Permutation: free
   Number of permutations: 999
   
   Model: capscale(formula = sqrt(scent.net.cut) ~ plants.majorpoll.cut)
            Df Variance      F Pr(>F)
   Model     3   808879 0.5422  0.792
   Residual 39 19393005

plot(scent.cap, type="text")

##Indicator compounds

library(indicspecies)

plants.majorpoll.cut <- plants.majorpoll[rownames(scent.net.cut)]
  
mp <- multipatt(as.data.frame(scent.net.cut), plants.majorpoll.cut, control=how(nperm=999))
summary(mp)

   
    Multilevel pattern analysis
    ---------------------------
   
    Association function: IndVal.g
    Significance level (alpha): 0.05
   
    Total number of species: 355
    Selected number of species: 11 
    Number of species associated to 1 group: 5 
    Number of species associated to 2 groups: 4 
    Number of species associated to 3 groups: 2 
   
    List of species associated to each combination: 
   
    Group butterfly  #sps.  2 
                                                                   stat p.value  
   1-Hexene, 3,5-dimethyl-                                        0.946   0.050 *
   Carbamic acid, N-phenyl-, 1,5-dimethyl-1-vinyl-4-hexenyl ester 0.942   0.038 *
   
    Group hummingbird  #sps.  3 
                                                   stat p.value  
   Propanoic acid, 2-methyl-, 3-methylbutyl ester 1.000   0.037 *
   Petasitene                                     0.996   0.037 *
   2H-Pyran-2-one, tetrahydro-3,6-dimethyl-       0.962   0.030 *
   
    Group butterfly+hummingbird  #sps.  3 
                                                              stat p.value   
   Bicyclo[7.2.0]undec-4-ene, 4,11,11-trimethyl-8-methylene- 0.967   0.002 **
   Methoxyacetic acid, tetradecyl ester                      0.941   0.015 * 
   1-Pentanol, 2-ethyl-4-methyl-                             0.918   0.050 * 
   
    Group fly+hummingbird  #sps.  1 
                         stat p.value  
   2-Butenal, 3-methyl- 0.921   0.043 *
   
    Group butterfly+fly+hummingbird  #sps.  2 
                          stat p.value  
   (E)-4-Oxohex-2-enal   0.941   0.012 *
   1-Propanol, 2-methyl- 0.890   0.039 *
   ---
   Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#Shannon diversity

###Do generalist plants have higher diversity of scent compounds?
plants.numlinks  <- rowSums(decostand(net.cut, "pa"))
#polls.classtree <- compute.brlen(polls.classtree$phylo,1)
#library(picante)
#plants.pdlinks   <- pd(net.lump, polls.classtree)$PD
plants.numcompounds <- rowSums(decostand(scent.net.cut, "pa"))
plants.sumcompounds <- rowSums(scent.net.cut)
plants.shannoncompounds <- diversity(scent.net.cut, index="shannon")
plants.shannonlinks <- diversity(net.grp, index="shannon")

scent.diversity <- tibble(numlinks=plants.numlinks, numcompounds=plants.numcompounds, sumcompounds=plants.sumcompounds, shannoncompounds=plants.shannoncompounds, shannonlinks=plants.shannonlinks)

ggplot(scent.diversity, aes(y=numcompounds, x=numlinks))+
    geom_point()+ geom_smooth()

   `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

anova(lm(numlinks~numcompounds, data=scent.diversity))

   Analysis of Variance Table
   
   Response: numlinks
                Df Sum Sq Mean Sq F value  Pr(>F)  
   numcompounds  1  792.3  792.30  3.8827 0.05556 .
   Residuals    41 8366.4  204.06                  
   ---
   Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#tells us if plants with more diverse scent get more pollinator species
anova(lm(plants.numlinks~shannoncompounds, data=scent.diversity))

   Analysis of Variance Table
   
   Response: plants.numlinks
                    Df Sum Sq Mean Sq F value Pr(>F)
   shannoncompounds  1   35.4  35.438  0.1593 0.6919
   Residuals        41 9123.3 222.520

#plot(plants.numcompounds~plants.pdlinks)

#plot(plants.shannoncompounds~plants.pdlinks)

shannon.compoundxlinks.plot <- ggplot(scent.diversity, aes(y=shannonlinks, x=shannoncompounds))+
    geom_point()+ geom_smooth()
ggsave(shannon.compoundxlinks.plot, file='shannon.compoundxlinks.plot.jpg', width=130, height=105, units= 'mm', dpi=300)

   `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

compoundxlinks.plot <- ggplot(scent.diversity, aes(y=plants.numlinks, x=shannoncompounds))+
    geom_point()+ geom_smooth()
ggsave(compoundxlinks.plot, file='compoundxlinks.plot.jpg', width=130, height=105, units= 'mm', dpi=300)

   `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

shannon.compound.plot <- ggplot(scent.diversity, aes(y=shannoncompounds, x=plants.majorpoll.cut))+
    geom_point()+ geom_smooth()
ggsave(shannon.compound.plot, file='shannon.compound.plot.jpg', width=110, height=105, units= 'mm', dpi=300)

   `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

numcompounds.plot <-ggplot(scent.diversity, aes(y=numcompounds, x=plants.majorpoll.cut))+
    geom_point()+ geom_smooth()
ggsave(numcompounds.plot, file='numcompounds.plot.jpg', width=110, height=105, units= 'mm', dpi=300)

   `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

sumcompounds.plot<- ggplot(scent.diversity, aes(y=sumcompounds, x=plants.majorpoll.cut))+
    geom_point()+ geom_smooth()
ggsave(sumcompounds.plot, file='sumcompounds.plot.jpg', width=110, height=105, units= 'mm', dpi=300)

   `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

#VOC-pollinator network

#interaction web plot
#voc-pollinator bipartite

#make dataframe with vocs and pollinators
##change scent.net to long format
plantcompoundemission <- scent.net.cut %>% as.data.frame() %>% rownames_to_column("plant") %>% #average peak area
  pivot_longer(-plant, names_to = "compound", values_to = "emissions") %>% filter(emissions!=0)


#write_csv(plantcompoundemission, "plantcompoundemissions.csv")
 
#total number of visits to plant by each pollinator species
plantpollinatorvisits <- net.cut %>% as.data.frame() %>% rownames_to_column("plant") %>% pivot_longer(-plant, names_to = "pollinator", values_to = "visits")


#total number of visits to each compound. Pollinators rows, compounds cols, values = number of visits
scentpollnet <- left_join(plantcompoundemission, plantpollinatorvisits,relationship = "many-to-many") %>%  group_by(pollinator, compound) %>% summarise(visits=sum(visits)) %>%  pivot_wider(names_from = "compound", values_from = "visits") %>% column_to_rownames("pollinator") %>% as.matrix()

   Joining with `by = join_by(plant)`
   `summarise()` has grouped output by 'pollinator'. You can override using the
   `.groups` argument.

#relative amount of each compound
plantcompoundrelative <- scent.net.cut %>% as.data.frame() %>% decostand(method="tot") %>% rownames_to_column("plant") %>% #average peak area
  pivot_longer(-plant, names_to = "compound", values_to = "emissions") %>% filter(emissions!=0)

#multiply number of visits by relative peak area
# weights visit to compounds. Proportions visits based on how many VOCs there are.
scentpollnet.area <- left_join(plantcompoundrelative, plantpollinatorvisits,relationship = "many-to-many") %>%  group_by(pollinator, compound) %>% mutate(visits.area=visits*emissions) %>%  summarise(visits.area=sum(visits.area)) %>%  pivot_wider(names_from = "compound", values_from = "visits.area") %>% column_to_rownames("pollinator") %>% as.matrix()

   Joining with `by = join_by(plant)`
   `summarise()` has grouped output by 'pollinator'. You can override using the
   `.groups` argument.

scentpollnet.class <- tibble(pollinator=rownames(scentpollnet)) %>% left_join(polls.class.dataframe) %>% pull(pollinator_group)

   Joining with `by = join_by(pollinator)`

#plot bipartite network
#plotweb(sqrt(scentpollnet), text.rot=90, method="cca", col.high=pal[2], bor.col.high=pal[2], #col.low=pcols[scentpollnet.class], bor.col.low=pcols[scentpollnet.class], col.interaction=pal[5], #bor.col.interaction=pal[5])

#plotweb for scentpollnet.area

plotweb(sqrt(scentpollnet.area), text.rot=90, method="cca", col.high=pal[2], bor.col.high=pal[2], col.low=pcols[scentpollnet.class], bor.col.low=pcols[scentpollnet.class], col.interaction=pal[5], bor.col.interaction=pal[5])

##Heatmap scent-net

#value = number of 
heatmaply(scentpollnet, scale="none", showticklabels = c(F,T), margins= c(0,NA,0,0),
          scale_fill_gradient_fun = ggplot2::scale_fill_gradient(low = "white", high = "#000050", trans="log1p"),
          distfun=vegdist, 
          hclustfun=function(x) hclust(x, method="mcquitty"), 
          #row_side_colors = data.frame(P=plants.majorpoll[rowSums(net)>0]),
          #row_side_palette=function(n) pcols, 
          row_side_colors=data.frame(PollinatorGroup=scentpollnet.class), 
          row_side_palette=function(n) pcols.all)

voc.pollgrp.cap <- capscale(sqrt(scentpollnet.area)~scentpollnet.class)
anova(voc.pollgrp.cap)

   Permutation test for capscale under reduced model
   Permutation: free
   Number of permutations: 999
   
   Model: capscale(formula = sqrt(scentpollnet.area) ~ scentpollnet.class)
            Df Variance    F Pr(>F)
   Model     4   15.904 1.67  0.105
   Residual 88  209.505

plot(voc.pollgrp.cap, type="text")

#chemical class x pollinator group bipartite

#make vector with compound class and all scent compounds
vol.class <-deframe(select(compound.class, compound_long, major_class) )

#turns scentpollnet matrix to data frame then pivots longer to have 2 cols called chemical and visits

chemical.taxa <- left_join(plantcompoundemission, plantpollinatorvisits,relationship = "many-to-many") %>% left_join(polls.class.dataframe) %>% left_join(compound.class, by=c("compound"= "compound_long")) %>% group_by(pollinator_group, major_class) %>% summarise(visits=sum(visits)) %>% pivot_wider(names_from = major_class, values_from = visits) %>% column_to_rownames("pollinator_group") %>% as.matrix()

   Joining with `by = join_by(plant)`
   Joining with `by = join_by(pollinator)`
   `summarise()` has grouped output by 'pollinator_group'. You can override using
   the `.groups` argument.

#chemical class & pollinator group network
# how do I incorporate chemical class instead of compounds? 
#make data frame with chemical class and pollinator family - use scent.net & chemical.class
#and use scentpollnet but scentpollnet compounds in different order than scent.net
plotweb(sqrt(chemical.taxa), text.rot=90, method="normal", col.high=pal[3], bor.col.high=pal[3], col.low=pcols.all, bor.col.low=pcols.all, col.interaction=pal[5], bor.col.interaction=pal[5])

library(bipartite)
library(vcd)

   Loading required package: grid

#side quest
mosaic(chemical.taxa)

mosaic(net)

chisq.test(chemical.taxa) #p=1 so this means that no specialist compound classes. and vise versa. nobody is specialized 

   
    Pearson's Chi-squared test
   
   data:  chemical.taxa
   X-squared = 8488.8, df = 24, p-value < 2.2e-16

#relative peak area x visits

#value = number of 
heatmaply(scentpollnet.area, scale="none", showticklabels = c(F,T), margins= c(0,NA,0,0),
          scale_fill_gradient_fun = ggplot2::scale_fill_gradient(low = "white", high = "#000050", trans="log1p"),
          distfun=vegdist, 
          hclustfun=function(x) hclust(x, method="mcquitty"), 
          #row_side_colors = data.frame(P=plants.majorpoll[rowSums(net)>0]),
          #row_side_palette=function(n) pcols, 
          row_side_colors=data.frame(PollinatorGroup=scentpollnet.class), 
          row_side_palette=function(n) pcols.all)

#Rarifaction

#species accumulation curves
#figure out what percent Pauls data catpured the whole pollinator community
#pull out 1 random 1 hr sampling period and say how many pollinators observered, then do this again and again randomly - variation in samples = standard error y axis= number of species, x= number of pollinator observations 

#class= bee, fly, butterfly, hummingbird
#xaxis is number of plant species observed (if you did the whole summer but only looked at 1 species and then next summer looked at 2 species, etc.)
#

#simplify net.long to only have numbers in it so that accumulation curve will work 
#rows are observation "units" and columns are pollinators
net.wide <- net.long %>% select(pollinator, interactions) %>% mutate(index=row_number()) %>% 
    pivot_wider(names_from=pollinator, values_from = interactions) %>% mutate(across(everything(),~replace_na(.x, 0))) %>% select(-index) %>% as.matrix()

for(class in unique(polls.class)) { #could replace class with family
  class_subset <- net[,polls.class==class] #making accumulation curve for each pollinator class
  class_sa <- vegan::specaccum(class_subset, permutations = 999) 
  plot.new()
  plot(class_sa, col = as.integer(factor(class))+1, add = class != "bee", lwd=2, 
       xlab="Pollinator Observations", ylab="Pollinator Species", main="Pollinator Species accumulation curve", sub=class)
}

legend("bottomright", unique(polls.class), 
       col = 2:5, pch = 19, title = "Class")

sa<- vegan::specaccum(net.wide,permutations=100, method="exact")
plot(sa) 

#Pauls total pollinator observations 8hrs/wk * 42wks = 336 hours total

visitationrate <- nrow(net.long)/336 #15.33 unique p-p observations per hr (observation can include multiple visits)
#observation=within a 15min period how many p-p interactions did we see, total number of links observed
2000/visitationrate # 130 hrs per treatment ( 2 treatments = 260 hrs per year)

   [1] 130.4854

#then based on my observation hours, what percent of the pollinator community could be captured

#normalized degree for p-p

library(rstatix)

   
   Attaching package: 'rstatix'

   The following object is masked from 'package:stats':
   
       filter

ND_poll<-ND(net, normalised = T)$higher 

#creates dataframe with normalized degree for each pollinator and pollinator class (same as data.frame)
ND.poll.class <- tibble(ND=ND_poll, polls.class=polls.class) 
  
kruskal_test(ND.poll.class, ND~polls.class)

   # A tibble: 1 × 6
     .y.       n statistic    df     p method        
   * <chr> <int>     <dbl> <int> <dbl> <chr>         
   1 ND       93      4.01     4 0.404 Kruskal-Wallis

ggplot(ND.poll.class, aes(x= polls.class, y=ND)) +geom_boxplot()

#normalized degree for VOC-pollinator

library(rstatix)

#ND for pollinators in scent network
ND_poll_scent<-ND(scentpollnet, normalised = T)$lower 

#ND of VOC chemical class
ND_scent_poll<-ND(scentpollnet, normalised = T)$higher


#creates dataframe with normalized degree for each pollinator and pollinator class (same as data.frame)
ND.poll <-  tibble(ND=ND_poll_scent, polls.class=polls.class)

  
  
kruskal_test(ND.poll, ND~polls.class)

   # A tibble: 1 × 6
     .y.       n statistic    df     p method        
   * <chr> <int>     <dbl> <int> <dbl> <chr>         
   1 ND       93      2.69     4 0.611 Kruskal-Wallis

ggplot(ND.poll, aes(x= polls.class, y=ND)) +geom_boxplot()

#creates dataframe with normalized degree for each VOC and chemical class (same as data.frame)
ND.scent <- tibble(ND=ND_scent_poll, vols.class=compound.class$major_class) 
#as.data.frame(list(ND = ND_scent_poll, vols.class = compound.class$major_class))
  
kruskal_test(ND.scent, ND~vols.class)

   # A tibble: 1 × 6
     .y.       n statistic    df     p method        
   * <chr> <int>     <dbl> <int> <dbl> <chr>         
   1 ND      354      7.83     6 0.251 Kruskal-Wallis

ggplot(ND.scent, aes(x= vols.class, y=ND)) +geom_boxplot()

#normalized visits - home cooked

#created a function similar to normalized degree(ND) but ND looks at presence/absence in the network where normalized visits (NV) looks at the normalized weighted visits/interactions. This accounts for variation in the amount each plant or compound is visited and give those with higher interaction frequency greater weight or contribution to the network  

normalized.visits <- function (web, normalised = TRUE) 
{

  websum <- sum(web) #takes sum of all interactions in the web
  dlower <- rowSums(web)
  dhigher <- colSums(web)
  Nlow <- Nhigh <- 2
  if (normalised) {
    Nlow <- websum 
    Nhigh <- websum
  }
  low <- dlower/Nlow #weighting the interactions by dividing rowsums by the sum of all links in the network
  names(low) <- rownames(web) #dividing colsums by the sum of all links in the network
  high <- dhigher/Nhigh
  names(high) <- colnames(web)
  list(lower = low, higher = high)
}

#weighted normalized visits for p-p (pollinators only)- home cooked

#Which pollinators have the most visits. of all the visits recorded, what percentage of them were by that pollinator (for a given plant what percentage of visits does it get out of the whole meadow, for a pollinator what percentage of visits does it give out of all visits)

#normalized visits for pollinators in scent net
NV_poll<-normalized.visits(net, normalised = T)$higher 

#creates dataframe with normalized degree for each pollinator and pollinator class (same as data.frame)
NV.poll.class <- tibble(NV=NV_poll, polls.class=polls.class) 
  
kruskal_test(NV.poll.class, NV~polls.class)

   # A tibble: 1 × 6
     .y.       n statistic    df      p method        
   * <chr> <int>     <dbl> <int>  <dbl> <chr>         
   1 NV       93      11.6     4 0.0208 Kruskal-Wallis

ggplot(NV.poll.class, aes(x= polls.class, y=NV)) +geom_boxplot()

#only a few species making most of the visits
#1 bee responsible for over 40% of visits

#normalized visits for scent network

#normalized visits for pollinators in scent net
#which pollinator class visits the most amount of compounds weighted by the number of visits
NV_poll_scentnet<-normalized.visits(scentpollnet, normalised = T)$lower 

#creates dataframe with normalized degree for each pollinator and pollinator class (same as data.frame)
NV.poll.class.scentnet <- tibble(NV=NV_poll_scentnet, polls.class=polls.class) 
  
kruskal_test(NV.poll.class.scentnet, NV~polls.class)

   # A tibble: 1 × 6
     .y.       n statistic    df     p method        
   * <chr> <int>     <dbl> <int> <dbl> <chr>         
   1 NV       93      2.91     4 0.573 Kruskal-Wallis

ggplot(NV.poll.class.scentnet, aes(x= polls.class, y=NV)) +geom_boxplot()

#super generalist bees

#normalized visits for VOCs in scent net
#which compound class are most visited
NV_voc<-normalized.visits(scentpollnet, normalised = T)$higher 

#creates dataframe with normalized degree for each pollinator and pollinator class (same as data.frame)
NV.vols.class <- tibble(NV=NV_voc, vols.class=compound.class$major_class) 
  
#TODO: Fix 
k_t.vols.class <- kruskal_test(NV.vols.class, NV~vols.class)


library(coin)

   Loading required package: survival

   
   Attaching package: 'coin'

   The following objects are masked from 'package:rstatix':
   
       chisq_test, friedman_test, kruskal_test, sign_test, wilcox_test

   The following object is masked from 'package:sna':
   
       rperm

library(PMCMRplus) #for posthoc
library(PMCMR)

   PMCMR is superseded by PMCMRplus and will be no longer maintained. You may wish to install PMCMRplus instead.

#posthoc_vols.class <- posthoc.kruskal.nemenyi.test(NV.vols.class$vols.class, NV.vols.class$NV, dist = "tukey")

ggplot(NV.vols.class, aes(x= vols.class, y=NV)) +geom_boxplot()

#tells you what vocs get the most visits
#no difference in visitation rate to compound classes

#modularity for net

library(vcd)

modules.net <- computeModules(sqrt(net))

plotModuleWeb(modules.net, labsize=.5)

module.net.list <- listModuleInformation(modules.net)

flatten(module.net.list[[2]])

   [[1]]
    [1] "androsace_septentrionalis" "claytonia_lanceolata"     
    [3] "draba_aurea"               "erythronium_grandiflorum" 
    [5] "hydrophyllum_capitatum"    "lomatium_dissectum"       
    [7] "pseudocymopterus_montanus" "ranunculus_inamoenus"     
    [9] "taraxacum_officinale"      "viola_praemorsa"          
   [11] "amelanchier_alnifolia"     "fragaria_virginiana"      
   [13] "hydrophyllum_fendleri"     "linum_lewisii"            
   
   [[2]]
    [1] "andrena_cyanophila"     "lasioglossum_spp"       "sphecodes_sp_01"       
    [4] "syrphidae_spp"          "scathophagidae_sp_01"   "halictus_confusus"     
    [7] "halictus_virgatellus"   "halictidae_green_spp_1" "eupeodes_volucris"     
   [10] "andrena_sp_02"          "muscidae_sp_02"         "syrphidae_02"          
   [13] "rhamphomyia_sp_01"      "cartosyrphus_tarda"     "callophrys_affinis"    
   [16] "osmia_lignaria"         "cartosyrphus_sp_02"     "tephritidae_01"        
   [19] "polygonia_zephyrus"    
   
   [[3]]
   [1] "delphinium_nuttallianum" "ipomopsis_aggregata"    
   [3] "gentiana_parryi"         "vicia_americana"        
   
   [[4]]
   [1] "bombus_appositus"        "selasphorus_platycercus"
   [3] "anthidium_emarginatum"   "bombus_nevadensis"      
   [5] "papilio_gothica"        
   
   [[5]]
    [1] "agoseris_glauca"      "arenaria_congesta"    "eriogonum_subalpinum"
    [4] "eriogonum_umbellatum" "galium_boreale"       "mahonia_repens"      
    [7] "potentilla_hippiana"  "rosa_woodsii"         "sedum_lanceolatum"   
   [10] "achillea_millefolium" "potentilla_gracilis" 
   
   [[6]]
    [1] "muscidae_spp_01"          "colletes_kincaidii"      
    [3] "halictus_rubicundus"      "halictidae_spp_01"       
    [5] "panurginus_ineptus"       "thricops_septentrionalis"
    [7] "osmia_bucephala"          "chrysis_coerulans"       
    [9] "villa_eumenes"            "sphaerophoria_robusta"   
   [11] "glaucopsyche_lygdamus"    "pieris_mcdunnoughi"      
   [13] "andrena_vicinoides"       "colletes_nigrifrons"     
   [15] "melanstoma_caerulescens" 
   
   [[7]]
    [1] "agoseris_aurantiaca"        "campanula_rotundifolia"    
    [3] "erigeron_speciosus"         "helianthella_quinquenervis"
    [5] "heliomeris_multiflora"      "heterotheca_villosa"       
    [7] "senecio_serra"              "solidago_multiradiata"     
    [9] "castilleja_sulphurea"       "phacelia_heterophylla"     
   
   [[8]]
    [1] "nymphalidae_spp"         "bombus_bifarius"        
    [3] "megachile_relativa"      "tachinidae_sp_01"       
    [5] "bombus_sylvicola"        "lycaenidae_spp"         
    [7] "anthophora_terminalis"   "bombus_flavifrons"      
    [9] "andrena_costillensis"    "bombus_californicus"    
   [11] "bombus_occidentalis"     "megachile_pugnata"      
   [13] "arctophila_flagrans"     "coelioxus_funeraria"    
   [15] "osmia_grindeliae"        "pieridae_spp"           
   [17] "speyeria_mormonia"       "chrysotoxum_ventricosum"
   [19] "megachile_melanophea"    "psithyrus_insularis"    
   [21] "lycaena_rubidus_sirius"  "colias_alexandra"       
   [23] "bombus_frigidus"         "mesebrina_latreillei"   
   [25] "syrphidae_spp_2"         "syrphidae_spp_4"        
   [27] "sphex_spp_1"             "bombus_balteatus"       
   [29] "syrphidae_spp_1"         "phormia_spp_1"          
   [31] "megachile_frigida"       "lycaena_cupreus"        
   [33] "bombus_mixtus"           "osmia_tristela"         
   [35] "melissodes_confusa"     
   
   [[9]]
   [1] "boechera_stricta"     "erigeron_flagellaris" "lathyrus_lanszwertii"
   [4] "mertensia_ciliata"    "mertensia_fusiformis" "senecio_integerrimus"
   
   [[10]]
    [1] "ochlodes_sylvanoides" "systoechus_sp"        "osmia_proxima"       
    [4] "pieris_rapa"          "andrena_transnigra"   "hoplitus_fulgida"    
    [7] "osmia_iridis"         "osmia_subaustralis"   "osmia_montana"       
   [10] "osmia_coloradensis"   "melanostoma_kelloggi" "syrphidae_spp_3"     
   [13] "eupeodes_lapponicus"  "euphydras_anicia"     "eristalis_latifrons" 
   [16] "osmia_sp_sgo"         "cryptopogon_sp_01"    "hoplitus_robusta"    
   [19] "coenonympha_ochracea"

pollinator.modules <- tibble(module=1:5) %>% mutate(pollinator=map(module, ~module.net.list[[2]][[.x]][[2]])) %>% unnest(pollinator) %>% left_join(polls.class.dataframe) %>% count(module, pollinator_group) %>% pivot_wider(names_from = pollinator_group, values_from = n, values_fill = 0) %>% column_to_rownames("module") %>% as.matrix()

   Joining with `by = join_by(pollinator)`

chisqr.polls <- chisq.test(pollinator.modules)

   Warning in chisq.test(pollinator.modules): Chi-squared approximation may be
   incorrect

residules.chisqr.polls <- residuals(chisqr.polls)

mosaic(t(pollinator.modules), shade = T, gp=shading_hcl)

#test to see p-p network modules
#mosaic(net, shade = T, gp=shading_hcl)

#modularity for scent net

modules.scentnet <- computeModules(sqrt(t(scentpollnet)))

#plotModuleWeb(modules.scentnet)

module.scentnet.list <- listModuleInformation(modules.scentnet)


pollinator.modules <- tibble(module=1:4) %>% mutate(pollinator=map(module, ~module.scentnet.list[[2]][[.x]][[2]])) %>% unnest(pollinator) %>% left_join(polls.class.dataframe) %>% count(module, pollinator_group) %>% pivot_wider(names_from = pollinator_group, values_from = n, values_fill = 0) %>% column_to_rownames("module") %>% as.matrix()

   Joining with `by = join_by(pollinator)`

chisqr.polls <- chisq.test(pollinator.modules)

   Warning in chisq.test(pollinator.modules): Chi-squared approximation may be
   incorrect

residules.chisqr.polls <- residuals(chisqr.polls)

mosaic(t(pollinator.modules), shade = T, gp=shading_hcl, labeling = vcd::labeling_border(rot_labels = c(0, 0)))

#chemical class modularity 


chemical.modules <- tibble(module=1:4) %>% mutate(compound_long=map(module, ~module.scentnet.list[[2]][[.x]][[1]])) %>% unnest(compound_long) %>% left_join(compound.class) %>% count(module, major_class) %>% pivot_wider(names_from = major_class, values_from = n, values_fill = 0) %>% column_to_rownames("module") %>% as.matrix()

   Joining with `by = join_by(compound_long)`

mosaic(t(chemical.modules), shade = T, gp=shading_hcl, labeling = vcd::labeling_border(rot_labels = c(0, 0)))

#stacked bar chart with flowers visited by each pollinator

library(ggplot2)
library(tidyverse)
library(RColorBrewer)
mypalette= c(brewer.pal(8, "Set2"), brewer.pal(12, "Set3"))

#looks at top 20 plants and top 75 pollinators (interactions sqrt to adjust scale on figure)
#pollinator perspective
ggplot(net.long %>% mutate(plant=fct_lump_n(plant, 18), pollinator=fct_lump_n(pollinator, 75)), aes(x=sqrt(interactions), y=pollinator, fill=plant)) + geom_col() + scale_fill_manual(values=mypalette)

#plant perspective
ggplot(net.long %>% mutate(plant=fct_lump_n(plant, 45), pollinator=fct_lump_n(pollinator, 18)), aes(x=sqrt(interactions), y=plant, fill=pollinator)) + geom_col() + scale_fill_manual(values=mypalette)

#phenology

phenology <- read_sheet("15Jnx8OWMzpakC0YvJWmZMmaRSnA0X-7HBjymDsvUfjY")

   ! Using an auto-discovered, cached token.

     To suppress this message, modify your code or options to clearly consent to
     the use of a cached token.

     See gargle's "Non-interactive auth" vignette for more details:

     <https://gargle.r-lib.org/articles/non-interactive-auth.html>

   ℹ The googlesheets4 package is using a cached token for 'jbohey@uci.edu'.

   ✔ Reading from "caradonna_rmbl_flowering_phenology_data_EDI".

   ✔ Range 'caradonna_rmbl_flowering_phenology_data_EDI'.

#stack plots and heatmaps 

#Flowers#
ggplot(phenology, aes(x=week_num, y=flower_count, fill= plant )) + geom_col() + scale_fill_manual(values=sample(rainbow(46))) +facet_wrap(vars(year))

#heatmap
ggplot(phenology, aes(x=week_num, y=plant, fill= sqrt(flower_count ))) + geom_tile() + scale_fill_viridis_c()  + facet_wrap(vars(year))

#pollinators#

ggplot(net.long, aes(x=week_num, y=interactions, fill= pollinator )) + geom_col() + scale_fill_manual(values=sample(rainbow(93))) +facet_wrap(vars(year))

#heatmap
ggplot(net.long, aes(x=week_num, y=pollinator, fill= sqrt(interactions ))) + geom_tile() + scale_fill_viridis_c()  + facet_wrap(vars(year))

#scent#
#number of flowers per week and multiple by mean scent of that species

#how many flowers are open each week
flowers.per.week <- phenology %>% group_by(year, week_num, plant) %>%  
                    summarise(flower_count=sum(flower_count))

   `summarise()` has grouped output by 'year', 'week_num'. You can override using
   the `.groups` argument.

#ideally want for each week, 1 average weighted mean of scent compound 
  #but we only have a few scent samples so just taking the mean of scent compounds for each plant

#flowers per week percentage 
flower.percentage <- flowers.per.week %>% group_by(year, week_num) %>% 
        mutate(flower_count=flower_count/sum(flower_count)) %>% 
        pivot_wider(names_from = c(year,week_num), values_from = flower_count, values_fill = 0) %>% 
        filter(plant%in%rownames(scent.net)) %>% 
        arrange(match(plant, rownames(scent.net)))%>%  #matches scent.net plant order
        column_to_rownames("plant")  %>% as.matrix()
#change names in pauls data

scentxphenology <- t(flower.percentage) %*% scent.net[rownames(flower.percentage),]  
#values in matrix = weighted average of flower scents and mean is weighted by how many flowers there are that week

#replaces Nans with 0
scentxphenology[is.nan(scentxphenology)] <- 0

#ggplot likes long format
scent.phenology <- scentxphenology %>% as.data.frame() %>%  rownames_to_column("yearweek") %>% separate(yearweek, into = c("year", "week_num")) %>% mutate(year=as.integer(year), week_num= as.integer(week_num)) %>%  pivot_longer(all_of(colnames(scent.net)))

ggplot(scent.phenology, aes(x=week_num, y=value, fill= name )) + geom_col() + scale_fill_manual(values=sample(rainbow(355)), guide="none") +facet_wrap(vars(year))

#heatmap
ggplot(scent.phenology, aes(x=week_num, y=name, fill= sqrt(value ))) + geom_tile() + scale_fill_viridis_c()  + facet_wrap(vars(year))

#mega table ##flower and poll traits ###VOCS

#number of compounds each plant species has
#total emissions for each species
#diversity of scents (shannon diversity)
library(tidyverse)

scent.net.long <- scent.net %>% as.data.frame() %>% rownames_to_column("plant") %>% #average peak area
  pivot_longer(-plant, names_to = "compound", values_to = "emissions") %>% filter(emissions!=0)


compound.class.percentage <- scent.net.long %>% left_join(compound.class, by=c(compound="compound_long")) %>%  group_by(plant, major_class) %>%  summarise(emissions=sum(emissions)) %>%  #gives total emissions for each compound class
        mutate(emissions=emissions/sum(emissions)) %>%  #relative amount of each compound class in species scent
        pivot_wider(names_from = major_class, values_from = emissions, names_prefix = "percent_", values_fill = 0) #makes column for each compound class and gives percentage of that compound class out of the total emissions for that species 

   `summarise()` has grouped output by 'plant'. You can override using the
   `.groups` argument.

#rare vs common compounds (mean rarity)
##for each compound, how many times does it occur in the dataset? 

compound.rarity <- scent.net.long %>%  count(compound) %>%  #number of species the compound occurs in
          mutate(rarity=n/nrow(scent.net))  #the percentage of species the compound occurs in
 
weighted.scent.rarity <- scent.net.long %>% left_join(compound.rarity) %>% group_by(plant) %>%    
            mutate(emissions=emissions/sum(emissions)) %>%  #makes emissions -> relative emissions    
            mutate(weighted_scent_rarity = emissions*rarity) %>% #multiplies emissions x rarity to get weighted average
            summarise(weighted_scent_rarity= sum(weighted_scent_rarity)) #then add the different parts of the average

   Joining with `by = join_by(compound)`

 #rarity score of 1 means this compound is found in all the plant species, 0= no other species has these compounds 
#mertensia has 0.70 rarity - the weighted average of all the compounds rarity; 70% of compounds (or total emissions) are shared, 30% of compounds or total emissions other species don't have



scent.summary <- tibble(plant=rownames(scent.net), total_emissions= rowSums(scent.net), number_compounds=rowSums(scent.net>0),scent_shannon_diversity= diversity(scent.net, index="shannon")) %>%
  left_join(weighted.scent.rarity) %>%     
  left_join(compound.class.percentage) 

   Joining with `by = join_by(plant)`
   Joining with `by = join_by(plant)`

##floral traits summary

#pollinators that visit plants

#Percentage of visits by each pollinator family to each plant species
pollinator.family.percentage <- net.long %>% group_by(plant, pollinator_family) %>%  summarise(interactions=sum(interactions)) %>%  #gives sum of interactions for each pollinator family
        mutate(interactions=interactions/sum(interactions)) %>%  #relative amount of each pollinator family visiting each flower species
        pivot_wider(names_from = pollinator_family, values_from = interactions, names_prefix = "percent_", values_fill = 0)

   `summarise()` has grouped output by 'plant'. You can override using the
   `.groups` argument.

#Percentage of visits by each pollinator group to each plant species
pollinator.group.percentage <- net.long %>% group_by(plant, pollinator_group) %>%  summarise(interactions=sum(interactions)) %>%  #gives sum of interactions for each pollinator group
        mutate(interactions=interactions/sum(interactions)) %>%  #relative amount of each pollinator group visiting each flower species
        pivot_wider(names_from = pollinator_group, values_from = interactions, names_prefix = "percent_", values_fill = 0)

   `summarise()` has grouped output by 'plant'. You can override using the
   `.groups` argument.

#rarity - how generalized the pollinator is
pollinator.rarity <- net.long %>% count(plant, pollinator) %>%  #adds how many rows of interactions between plant & poll species
          count(pollinator, name="number_plants_visited") %>%  #number of plants the pollinator goes to
          mutate(rarity=number_plants_visited/nrow(net))  #the percentage of plant species the pollinator visits in the meadow

#gives for each plant species, which pollinators visited and their relative contributions to total number of visits
 net.long.relative.interactions <- net.long %>% group_by(plant, pollinator) %>% 
            summarise(interactions=sum(interactions)) %>% 
            left_join(pollinator.rarity) %>% 
            group_by(plant) %>%    
            mutate(interactions=interactions/sum(interactions))

   `summarise()` has grouped output by 'plant'. You can override using the
   `.groups` argument.
   Joining with `by = join_by(pollinator)`

 #percent of interactions that come from each pollinator species    
 weighted.pollinator.rarity <-   net.long.relative.interactions %>% mutate(weighted_poll_rarity = interactions*rarity) %>% #multiplies interactions x rarity to get weighted average
            summarise(weighted_poll_rarity = sum(weighted_poll_rarity)) 

flower.visitor.summary <- tibble(plant=rownames(net), total_visits= rowSums(net), number_pollspecies_visited=rowSums(net>0),pollinator_shannon_diversity= diversity(net, index="shannon")) %>%
  left_join(weighted.pollinator.rarity) %>% 
  left_join(pollinator.family.percentage) %>% left_join(pollinator.group.percentage)

   Joining with `by = join_by(plant)`
   Joining with `by = join_by(plant)`
   Joining with `by = join_by(plant)`

megatable <- scent.summary %>% full_join(flower.visitor.summary) #full join bc pauls data and scent data plants don't match

   Joining with `by = join_by(plant)`

#correlations

library(pheatmap)
library(smatr)


correlation.table <- megatable %>% column_to_rownames("plant") %>% cor(use="na.or.complete")

pheatmap(correlation.table, scale="none", clustering_method = "mcquitty", color=colorRampPalette(c("blue","white", "red"))(200),
         breaks = seq(-1,1, by=0.01)) #which ones it labels on the scale

#scatterplot with just weighted scent and weighted pollinator visits
weighted.scent.poll <- ggplot(megatable, aes(x=weighted_poll_rarity, y=weighted_scent_rarity, color=number_pollspecies_visited)) + geom_point()+geom_smooth() + geom_text(aes(label=plant))     

print(weighted.scent.poll)

   `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

   Warning: Removed 5 rows containing non-finite outside the scale range
   (`stat_smooth()`).

   Warning: The following aesthetics were dropped during statistical transformation:
   colour.
   ℹ This can happen when ggplot fails to infer the correct grouping structure in
     the data.
   ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
     variable into a factor?

   Warning: Removed 5 rows containing missing values or values outside the scale range
   (`geom_point()`).

   Warning: Removed 5 rows containing missing values or values outside the scale range
   (`geom_text()`).

#the most common volatiles get the most visits from the most generalized pollinators

sma.weighted.scent.poll <- sma(weighted_poll_rarity~weighted_scent_rarity, data=megatable )

summary(sma.weighted.scent.poll)

   Call: sma(formula = weighted_poll_rarity ~ weighted_scent_rarity, data = megatable) 
   
   Fit using Standardized Major Axis 
   
   ------------------------------------------------------------
   Coefficients:
                elevation    slope
   estimate    -0.3602556 1.352232
   lower limit -0.5885207 1.001742
   upper limit -0.1319906 1.825352
   
   H0 : variables uncorrelated
   R-squared : 0.06769348 
   P-value : 0.091989

#percent_hesperiidae x nitrogenous vocs
ggplot(megatable, aes(x=percent_nitrogenous, y=percent_lycaenidae, color=number_pollspecies_visited)) + geom_point()+geom_smooth() + geom_text(aes(label=plant))

   `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

   Warning: Removed 5 rows containing non-finite outside the scale range
   (`stat_smooth()`).

   Warning: The following aesthetics were dropped during statistical transformation:
   colour.
   ℹ This can happen when ggplot fails to infer the correct grouping structure in
     the data.
   ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
     variable into a factor?

   Warning: Removed 5 rows containing missing values or values outside the scale range
   (`geom_point()`).

   Warning: Removed 5 rows containing missing values or values outside the scale range
   (`geom_text()`).

##pollinator trait summary